function [Param, Covar, Resid, Info ] = mvnrmle(Data, Design, ...
   MaxIter, TolParam, TolObj, Covar0, CovarFormat)
%MVNRMLE - Multivariate normal regression (ignore missing data).
% Estimate a NUMPARAMS x 1 column vector of model parameters Param and a NUMSERIES x NUMSERIES
% matrix of covariance parameters Covar for a multivariate normal regression model without missing
% data, where the model has the form
%		Data(k)  ~  N(Design(k) * Param, Covar)
% for samples k = 1, ... , NUMSAMPLES. If no output arguments, display a plot of convergence (or
% non-convergence) of the algorithm to a maximum likelihood estimate.
%
%		[Param, Covar] = mvnrmle(Data, Design);
%		[Param, Covar, Resid, Info ] = mvnrmle(Data, Design, ...
%			MaxIter, TolObj, TolParam, Covar0, CovarFormat);
%
%		mvnrmle(Data, Design);
%		mvnrmle(Data, Design, MaxIter, TolObj, TolParam, Covar0, CovarFormat);
%
% Inputs:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a NUMSERIES-dimensional random
%		vector. If a data sample has missing values (with NaNs), the sample is ignored (use
%		ECMMVNRMLE to handle missing data).
%	Design - Either a matrix or a cell-array to handle two distinct model structures. First, if
%		NUMSERIES = 1, Design can be a NUMSAMPLES x NUMPARAMS matrix with known values. This is the
%		"standard" form for regression on a single data series. Alternatively, for any
%		NUMSERIES >= 1, Design can be a cell array of either 1 or NUMSAMPLES cells, where each cell
%		contains a NUMSERIES x NUMPARAMS matrix of known values. If Design has a single cell, then
%		it is assumed to be the same Design matrix for each sample. Otherwise, Design must contain
%		individual Design matrices for each sample.
%
% Optional Inputs:
%	MaxIter - Maximum number of iterations for the estimation algorithm. Default value is 100.
%	TolParam - Convergence tolerance for estimation algorithm based on changes in model parameter
%		estimates. Default value is sqrt(eps) which is about 1.0e-8 for double precision. The
%		convergence test for changes in model parameters is
%			||Param(k) - Param(k-1)|| < TolParam * (1 + ||Param(k)||)
%		for iteration k = 2, 3, ... . Convergence is assumed when both the TolParam and TolObj
%		conditions are satisfied. If both TolParam <= 0 and TolObj <= 0, then do maximum number of
%		iterations (MaxIter) regardless of the results of convergence tests.
%	TolObj - Convergence tolerance for estimation algorithm based on changes in the objective
%		function. Default value is eps^(3/4) which is about 1.0e-12 for double precision. The
%		convergence test for changes in the objective function is
%			|Obj(k) - Obj(k-1)| < TolObj * (1 + |Obj(k)|)
%		for iteration k = 2, 3, ... . Convergence is assumed when both the TolParam and TolObj
%		conditions are satisfied. If both TolParam <= 0 and TolObj <= 0, then do maximum number of
%		iterations (MaxIter) regardless of the results of convergence tests.
%	Covar0 - NUMSERIES x NUMSERIES matrix that contains a user-supplied initial or known estimate
%		for the covariance matrix of the residuals of the regression. Default is an identity matrix.
%	CovarFormat - String that specifies the format for the covariance matrix. The choices are:
%		'full' - Default method. Compute the full covariance matrix.
%		'diagonal' - Force the covariance matrix to be a diagonal matrix.
%
% Outputs:
%	Param - NUMPARAMS x 1 column vector of estimates for the model parameters of the regression
%		model.
%	Covar - NUMSERIES x NUMSERIES matrix of estimates for the covariance of the residuals of the
%		regression.
%	Resid - NUMSAMPLES x NUMSERIES matrix of residuals from the regression. For any row with missing
%		values in Data, the corresponding row of residuals is represented as all NaN missing values
%		since this routine ignores rows with NaN values.
%	Info - A structure that contains additional information from the regression. The structure has
%		the following fields:
%		Info.Obj - A variable-extent column vector with no more than MAXITER elements that contains
%			each value of the objective function at each iteration of the estimation algorithm. The
%			last value in this vector, which is Obj(end), is the terminal estimate of the objective
%			function. If doing maximum likelihood estimation, the objective function is the
%			log-likelihood function.
%		Info.PrevParam - NUMPARAMS x 1 column vector of estimates for the model parameters from the
%			iteration just prior to the terminal iteration.
%		Info.PrevCovar - NUMSERIES x NUMSERIES matrix of estimates for the covariance parameters
%			from the iteration just prior to the terminal iteration.
%
% Notes:
%	This function does not accept an initial parameter vector since the parameters are estimated
%	directly from the first iteration onward.
%
%	Although the Design array should not have NaN values, ignored samples due to NaN values in Data
%	are also ignored in the corresponding Design array.
%
%	Note, however, that no NaN values are permitted in Design if it is a 1 x 1 cell array to be used
%	as a single Design matrix for each sample. In addition, a model in this form must have
%	NUMSERIES >= NUMPARAMS and rank(Design{1}) = NUMPARAMS.
%
%	Finally, note that the missing data routines ecmmvnrmle and ecmlsrmle are stricter about the
%	presence of NaN values in the Design array.
%
%	The optional output structure Info with fields Obj, PrevParam, and PrevCovar, is available for
%	diagnostic purposes.
%
% References:
%	[1] Roderick J. A. Little and Donald B. Rubin, Statistical Analysis with Missing Data, 2nd ed.,
%		John Wiley & Sons, Inc., 2002.
%	[2] Xiao-Li Meng and Donald B. Rubin, "Maximum Likelihood Estimation via the ECM Algorithm,"
%		Biometrika, Vol. 80, No. 2, 1993, pp. 267-278.
%
%	See also MVNRSTD, MVNROBJ, ECMMVNRMLE.

%	Copyright 2005-2007 The MathWorks, Inc.
%	$Revision: 1.1.6.4 $ $Date: 2007/07/03 20:42:11 $

% Step 1 - check arguments

if nargin < 7 || isempty(CovarFormat)
   CovarFormat = 'full';
else
   if ~any(strcmpi(CovarFormat,{'full','diagonal'}))
      error('Finance:mvnrmle:InvalidCovarianceFormat', ...
         'Invalid format specified for covariance matrix.');
   end
end
if nargin < 6
   Covar0 = [];
end
if nargin < 5 || isempty(TolObj)
   TolObj = eps^(3/4);
end
if nargin < 4 || isempty(TolParam)
   TolParam = sqrt(eps);
end
if nargin < 3 || isempty(MaxIter)
   MaxIter = 100;
elseif MaxIter < 1
   MaxIter = 1;
end
if nargin < 2
   error('Finance:mvnrmle:MissingInputArg', ...
      'Missing required arguments Data or Design.');
end

[NumSamples, NumSeries, NumParams] = checkmvnrsetup(Data, Design, [], Covar0, true);

if iscell(Design) && (numel(Design) == 1)
   SingleDesign = true;
else
   SingleDesign = false;
end

% Step 2 - observability and ignorability tests

Count = sum(any(isnan(Data),2));
if ((NumSamples - Count) * NumSeries) <= max(NumParams, (NumSeries * (NumSeries + 1))/2)
   if ((NumSamples - Count) * NumSeries) <= NumParams
      error('Finance:mvnrmle:InsufficientData', ...
         'Insufficient data to estimate either full or least-squares models.');
   else
      warning('Finance:mvnrmle:TryingReducedModel', ...
         'Insufficient data to estimate full model. Will try least-squares.');
      MaxIter = 1;
   end
end

% Step 3 - initialization

Param = [];

if isempty(Covar0)
   Covar = eye(NumSeries, NumSeries);
else
   Covar = Covar0;
end

if strcmpi(CovarFormat,'diagonal');		% need to make sure initial covar is diagonal if
   Covar = diag(diag(Covar));			% CovarFormat = 'diagonal' or else premature mle
end

Obj = [];

[CholCovar, CholState] = chol(Covar);
if CholState > 0
   warning('Finance:mvnrmle:NonPosDefCov', ...
      'Initial covariance not positive-definite. Will use identity matrix.');
   Covar = eye(NumSeries, NumSeries);
end

Resid = nan(NumSamples, NumSeries);

% Step 4 - main loop

for Iter = 1:MaxIter

   if nargout > 3
      Info.PrevParam = Param;
      Info.PrevCovar = Covar;
   end

   [CholCovar, CholState] = chol(Covar);

   LeftSide = zeros(NumParams, NumParams);
   RightSide = zeros(NumParams, 1);

   % Step 5 - parameter estimation step

   if iscell(Design)
      if SingleDesign
         A = CholCovar' \ Design{1};
         for k = 1:NumSamples
            % LeftSide kept in loop to maintain equivalence between
            % alternative SingleDesign forms
            if ~any(isnan(Data(k,:)))
               B = CholCovar' \ Data(k,:)';
               LeftSide = LeftSide + (A' * A);
               RightSide = RightSide + (A' * B);
            end
         end
      else
         for k = 1:NumSamples
            if ~any(isnan(Data(k,:)))
               A = CholCovar' \ Design{k};
               B = CholCovar' \ Data(k,:)';
               LeftSide = LeftSide + (A' * A);
               RightSide = RightSide + (A' * B);
            end
         end
      end
   else
      for k = 1:NumSamples
         if ~isnan(Data(k))
            A = CholCovar' \ Design(k,:);
            B = CholCovar' \ Data(k);
            LeftSide = LeftSide + (A' * A);
            RightSide = RightSide + (A' * B);
         end
      end
   end

   Param = LeftSide \ RightSide;

   % Step 6 - covariance estimation step

   Count = 0;
   Covar = zeros(NumSeries, NumSeries);
   if iscell(Design)
      if SingleDesign
         B = Design{1} * Param;
         for k = 1:NumSamples
            if ~any(isnan(Data(k,:)))
               Count = Count + 1;
               Resid(k,:) = Data(k,:) - B';
               Covar = Covar + (Resid(k,:)' * Resid(k,:));
            end
         end
      else
         for k = 1:NumSamples
            if ~any(isnan(Data(k,:)))
               Count = Count + 1;
               Resid(k,:) = Data(k,:) - (Design{k} * Param)';
               Covar = Covar + (Resid(k,:)' * Resid(k,:));
            end
         end
      end
   else
      for k = 1:NumSamples
         if ~isnan(Data(k))
            Count = Count + 1;
            Resid(k) = Data(k) - (Design(k,:) * Param);
            Covar = Covar + Resid(k)^2;
         end
      end
   end

   if strcmpi(CovarFormat,'diagonal')
      Covar = diag(diag(Covar));
   end

   Covar = (1.0/Count) .* Covar;

   % Step 7 - estimate objective and do convergence test

   Objective = mvnrobj(Data, Design, Param, Covar, CovarFormat);
   Obj = [Obj; Objective];

   if Iter > 1
      TestObj = Objective - Objective0;
      TestParam = norm(Param - Param0)/sqrt(NumParams);

      EpsObj = TolObj * (1 + abs(Objective));
      EpsParam = TolParam * (1 + norm(Param));

      if ((TestObj >= 0.0) && (TestObj < EpsObj)) && (TestParam < EpsParam)
         break
      end
   end

   Objective0 = Objective;
   Param0 = Param;

   if (Iter == MaxIter) && (MaxIter > 1) && (TolObj > 0.0) && (TolParam > 0.0)
      warning('Finance:mvnrmle:EarlyTermination', ...
         'Maximum iterations completed. Convergence criterion not satisfied.');
   end

   if NumSeries == 1
      break
   end
end

if nargout > 3
   Info.Obj = Obj;
end

% Step 8 - plot objective function if no output args

if ~nargout
   % Search for existing figure
   findfig = findall(0, 'Type', 'figure', 'tag', 'mvnrmle');

   if ~isempty(findfig)
      % Use existing figure
      f = figure(findfig);

   else
      % Create a new figure
      f = figure('Visible', 'Off', 'Tag', 'mvnrmle');
   end

   % Plot into current figure
   plot(Obj,'-bo','MarkerFaceColor','b','MarkerSize',3);
   title('\bfProgress of Estimation Algorithm in mvnrmle');
   xlabel('\bfIteration');
   ylabel('\bfLog-Likelihood');

   % Turn on visibility (if it was off)
   set(f, 'Visible', 'On')
end
